home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Amiga Collections: Topik
/
Topik - Disk 39 - Educational (19xx)(Topik Public Domain)(PD)[WB].zip
/
Topik - Disk 39 - Educational (19xx)(Topik Public Domain)(PD)[WB].adf
/
Planets
/
planet.c
< prev
next >
Wrap
C/C++ Source or Header
|
1991-03-09
|
16KB
|
664 lines
/* Planets
* A program to determine the position of
* the Sun., Mer., Ven., Mars, Jup, Sat & Ura
*
* reference Jean Meeus's book, " Astronomical Formulae for
* Calculators
* program by
*
* F.T. Mendenhall
* ihnp4!inuxe!fred
*
* Copyright 1985 by F.T. Mendenhall
* all rights reserved
*
* (This copyright prevents you from selling the program - the
* author grants anyone the right to copy and install the program
* on any machine it will run on)
!
! Modified by Alan Paeth, awpaeth@watcgl, January 1987 for
! (1) non-interactive mode (uses current GMT)
! (2) output in simplified Yale format for use with star charter software
! (3) lines 365-380 rewritten as a simplier C expression (awpaeth, Apr '87)
!
! Modified by Petri Launiainen, pl@intrin.FI, February 1987 for
! easier LOGFILE definition/modification in Makefile
!
! Modified by Alan Paeth, September 1987, for command line flags.
! The program continues to provide reduced Yale output in the old style,
! but new versions of "starchart" accept either. The old format allows
! representation of magnitudes below -.99 (useful for planets), but no
! spectral or stellar identification, which is not useful for planets.
!
! Modified by Jim Cobb, March 1988, to compute moon positions as well.
! Also, to change slightly the computation of epli, the obliquity
! of the ecliptic. Formerly epli had a term
! +0.00256*cos(omeg*rad).
! This term is a compensation for computing the apparent position
! of the sun (Chapter 18 in Meeus). But for converting
! ecliptical coordinates to right ascension and declination we want
! the mean value (that is, we don't want this term).
!
! Modified by Bob Leivian, Sept 1989, to compute local azimuth and altitude
! given a lat/lon, also fixed time so you can replace a single item
! without having to supply all (i.e. for today at 9:00 put -t 9)
! also ran the code through 'cb' to standardize the indention.
! I split the code in to more managible pieces.
! Also printed out time in local and UTC as well as JD
!
! See "Time Notes" section below if you plan to rewrite the user interface.
*/
#include <stdio.h>
#include <math.h>
/* default latitude longitude (of Tempe, AZ) -- replace with your favorite place */
#define DEFLONG -111.94
#define DEFLAT 33.32
#define DEFZONE 7 /* MST */
#ifndef PLANETFILE
#define PLANETFILE "./planet.star"
#endif
double pie, rad;
double htod(), atof(), kepler(), truean();
double longi(), lati(), poly(), to_int(), range();
int cls(), helio_trans(), geo_trans(), speak();
double planet_pos();
double sidereal_time();
int zone;
double longitude = DEFLONG;
double latitude = DEFLAT;
double GHA_Aries;
FILE *logfile;
char *progname;
/*
* Times Notes.
*
* This software would ideally use a generic UNIX system call to which would
* provide day and date info as integers, plus current time west of GMT.
* These values would become the default settings for each of the -m -d -y -t
* & -z. Unfortunately, time on across many UNIX systems does not (to my
* knowledge) exist as a uniform subroutine call.
*
* Instead, we presume that the very low level "gettimeofday" (GMT in seconds,
* from 1 Jan 1970) exists on all UNIX installations, and use this subroutine
* to use the present time of day when no command line parameters are
* present, as defaults for the current time have not been filled in.
*
* When any flag-style command line parameters are present, hardcoded defaults
* are substitude for -m, -d, -y -t switches, which may then be overridden.
* The -z flag is filled in by a companion return value from "gettimeofday".
* This approach makes for poor (on account of hard coding) defaults.
*
* A worthwhile change would be to simplify the command line control structure
* and default mechanism by using a higher level system call, or (safer),
* write the routine to convert Unix GMT into day, month, year, so that it
* lives within this program. This would allow for more reasonable defaults,
* and still merely one UNIX system call to get GMT time. Perhaps sources
* from a UNIX system may be studied to do this correctly.
*
* Expertise in "time" on a specific UNIX system is *NOT* what is required --
* what *IS* required is knowledge and access to enough independent UNIX
* systems to insure (to a high probability) that the code remains portable.
* In other words, only low-level system calls known to work across a large
* range of systems should be employed.
*
*/
#include <sys/time.h> /* for getting current GMT */
double
current_time(m, d, y, t, z)
int *m;
int *d;
int *y;
double *t;
int *z;
{
#define SECSPERDAY (60.0 * 60.0 * 24.0)
#ifdef UNIX
#define GMT1970 2440587.5
struct timeval tv;
struct timezone tz;
struct tm *p;
gettimeofday(&tv, &tz);
p = localtime(&tv);
*y = p->tm_year + 1900; /* local year */
*m = p->tm_mon + 1; /* local month */
*d = p->tm_mday; /* local day of month */
*t = p->tm_hour + (p->tm_min / 60.0); /* local time (0.0-23.999) */
*z = tz.tz_minuteswest; /* local time hrs west of Greenwich */
return (GMT1970 + tv.tv_sec / SECSPERDAY);
#endif
#ifdef AMIGA
time_t amiga_now;
struct tm *p;
#define GMT1978 2443509.5
amiga_now = time(&amiga_now);
p = localtime(&amiga_now);
*t = p->tm_hour + (p->tm_min / 60.0); /* local time (0.0-23.999) */
*z = DEFZONE * 60;
*y = p->tm_year + 1900; /* local year */
*m = p->tm_mon + 1; /* local month */
*d = p->tm_mday; /* local day of month */
return GMT1978 + (amiga_now / SECSPERDAY) + (DEFZONE / 24.0);
#endif
#ifdef NO_TIME_AVAIL
#define MST1989 2447528.291667
/* no time is avail, pick a default i.e. Noon Jan 1 1989 */
*y = 1989;
*m = 1;
*d = 1;
*t = 12.0;
*z = DEFZONE * 60.;
return MST1989;
#endif
}
double
commandtime(argc, argv)
char **argv;
{
int mm, dd, yy, aa1, bb1;
double jd, year, month, day, time, timez;
int j;
double now;
static char *usage =
"\nusage:planet {current time used}\nor\tplanet juliandate\nor\tplanet [-t hh.mm -z zone -y year -m mon -d day ]";
/* get current time for defaults */
now = current_time(&mm, &dd, &yy, &time, &zone);
timez = zone / 60.; /* time zone in hours */
/* no args - use current time */
if (argc == 1)
return now;
if ((argv[1][0] != '-')) {
/* one numeric arg - use it as the julian date */
if (argc == 2)
return (atof(argv[1]));
die("no switches - %s", usage);
} else {
/* flags present, set up defaults */
/* process any user overrides */
for (j = 1; j < argc; j++) {
if (argv[j][0] != '-')
die("unknown argument - %s", usage /*argv[j]*/);
switch (argv[j][1]) {
case 't':
time = htod(argv[++j]);
break;
case 'z':
timez = htod(argv[++j]);
if (timez < 0 || timez > 24)
die("time zone must be 0 to 24 west of UTC\n");
zone = timez * 60.;
break;
case 'y':
yy = atoi(argv[++j]);
break;
case 'm':
mm = atoi(argv[++j]);
if (mm < 1 || mm > 12)
die("month must be 1..12\n");
break;
case 'd':
dd = atoi(argv[++j]);
if (dd < 1 || dd > 31)
die("day must be 1..31\n");
break;
case 'l':
switch (argv[j][2]) {
case 'o':
longitude = atof(argv[++j]);
if (longitude < -180 || longitude > 180)
die("longitude must be +-180\n");
break;
case 'a':
latitude = atof(argv[++j]);
if (latitude < -90 || latitude > 90)
die("latitude must be +-90\n");
break;
default:
die("specify lat or lon\n");
}
break;
default:
die("unknown switch - %s", usage /*argv[j]*/ );
}
if (j == argc)
die("trailing command line flag - %s", argv[j - 1]);
}
}
day = ((float) dd + (time + timez) / 24.0);
if (mm > 2) {
year = yy;
month = mm;
} else {
year = yy - 1;
month = mm + 12;
}
aa1 = (int) (year / 100);
bb1 = 2 - aa1 + (int) (aa1 / 4);
jd = to_int(365.25 * year) + to_int(30.6001 * (month + 1));
jd = jd + day + 1720994.5;
if ((year + month / 100) > 1582.10)
jd = jd + bb1;
return (jd);
}
caltim(jd_frac, hour, min)
double jd_frac;
int *hour;
int *min;
{
double x;
x = jd_frac * 24.0;
*hour = (int) floor(x);
x -= floor(x);
x *= 60.0;
*min = (int) floor(x);
}
main(argc, argv)
int argc;
char **argv;
{
double jd, T0, epli;
double temp;
long int_jd;
int y, m, d, h, min, lh, lmin;
char *str;
progname = argv[0];
#ifdef AMIGA
logfile = 0; /* don't write a log */
#else
if ((logfile = fopen(PLANETFILE, "w")) == 0)
/*if ((logfile = fopen("/tmp/planet.star", "w")) == 0) {
die("fail on opening logfile %s\n", "/tmp/planet.star");
}*/
#endif
cls();
pie = 3.1415926535898;
rad = pie / 180.0;
/* calculate time T0 from 1900 January 0.5 ET */
jd = commandtime(argc, argv);
/* jd to calendar day */
int_jd = jd +.5001;
caldat(int_jd, &m, &d, &y);
/* julian day to UTC */
temp = jd - .4999;
caltim(temp - floor(temp), &h, &min);
printf("At %d/%d %d %d:%02d UTC (Julian: %.3f)\n", m, d, y, h, min, jd);
/* UTC to Local */
temp -= (zone / 1440.);
caltim(temp - floor(temp), &lh, &lmin);
int_jd = jd + .5001 - (zone / 1440.);
caldat(int_jd, &m, &d, &y);
if (lh > 12) {
lh -= 12;
str = "pm";
} else {
str = "am";
if (lh == 0)
lh = 12;
}
printf("local %d/%d %d:%02d%s (", m, d, lh, lmin, str);
/* compute exact local time for a given longitude */
temp = jd - .5 + (longitude / 360.);
caltim(temp - floor(temp), &lh, &lmin);
GHA_Aries = 15. * sidereal_time(jd);
printf("%2d:%02d true local at lat:%1.1f lon:%1.1f, GHA Aries:%1.1f)\n",
lh, lmin, latitude, longitude, GHA_Aries);
T0 = (jd - 2415020.0) / 36525.0;
printf("\nOBJECT R.A. DEC DISTANCE altitude azimuth\n");
epli = planet_pos(jd, T0);
moon_pos(T0, epli);
putchar('\n');
if (logfile)
close(logfile);
logfile = 0;
/* show local position of some stars also */
speak(101.1714, -16.7013, 0., -1.46, " " , "Sirius");
speak(213.7956, 19.236, 0., -0.04, " " , "Arcturus");
speak(88.6508, 7.4057, 0., 0.5, " " , "Betelgeuse");
exit(0);
} /* end of program main */
int
cls()
{
#define CLEARCNT 1
int i;
for (i = 0; i < CLEARCNT; i++)
putchar('\n');
}
double
poly(a0, a1, a2, a3, t)
double a0, a1, a2, a3, t;
{
return (a0 + a1 * t + a2 * t * t + a3 * t * t * t);
}
double
to_int(z)
double z;
{
long trunk;
trunk = (long) z;
z = (double) trunk;
return (z);
}
double
kepler(e, M)
double e, M;
{
double corr, e0, E0, E1;
e0 = e / rad;
corr = 1;
E0 = M;
while (corr > 0.000001) {
corr = (M + e0 * sin(E0 * rad) - E0) / (1 - e * cos(E0 * rad));
E1 = E0 + corr;
if (corr < 0) {
corr = -1.0 * corr;
}
E0 = E1;
}
return (E1);
}
double
truean(e, E)
double e, E;
{
double nu;
nu = 2.0 * atan(sqrt((1 + e) / (1 - e)) * tan(E * rad / 2));
nu = nu / rad;
if (nu < 0.0) {
nu = nu + 360.0;
}
if (fabs(nu - E) > 90.0) {
if (nu > 180.0) {
nu = nu - 180.0;
} else {
nu = nu + 180.0;
}
}
return (nu);
}
double
longi(w2, i, u)
double w2, i, u;
{
double x, y, l;
y = cos(i * rad) * sin(u * rad);
x = cos(u * rad);
l = atan2(y, x);
l = l / rad;
if (l < 0.0) {
l = l + 360.0;
}
return (l + w2);
}
double
lati(u, i)
double u, i;
{
double b;
b = asin(sin(u * rad) * sin(i * rad)) / rad;
return (b);
}
int
speak(ra, dec, dis, mag, sym, str)
double ra, dec, dis, mag;
char *sym, *str;
{
int h, m, s, deg;
double lha, sha, altitude, azimuth;
;
if (ra < 0) {
ra = ra + 360.0;
}
sha = 360. - ra;
ra = ra / 15.0; /* convert from degs to hours */
h = (int) to_int(ra);
m = (int) ((ra - h) * 60);
s = (int) (((ra - h) * 60 - m) * 60);
printf("%-12s%2dh %2dm %2ds ", str, h, m, s);
if (logfile)
fprintf(logfile, "%02d%02d%02d", h, m, s);
deg = (int) to_int(dec);
m = (int) ((dec - deg) * 60);
s = (int) (((dec - deg) * 60 - m) * 60);
if (m < 0) {
m = m * -1;
}
if (s < 0) {
s = s * -1;
}
#ifdef AMIGA
/* its printf doesn't know about signed */
if (deg > 0)
printf(" +%2ddeg %2dm %2ds ", deg, m, s);
else
printf(" %3ddeg %2dm %2ds ", deg, m, s);
#else
printf(" %+3ddeg %2dm %2ds ", deg, m, s);
#endif
if (dis > 0.0001)
printf(" %6.3f", dis);
else
printf(" ");
if (logfile)
if (mag < 0.0)
fprintf(logfile, "%+03d%02d-%2d%s %s\n",
deg, m, -(int) (10 * mag), sym, str);
else
fprintf(logfile, "%+03d%02d%3d%s %s\n",
deg, m, (int) (100 * mag), sym, str);
/* if we have a valid angle */
if (GHA_Aries > 0) {
/* now tell where it is at a given place on earth */
lha = range(GHA_Aries + sha + longitude);
altitude = sin(latitude * rad) * sin(dec * rad) +
cos(latitude * rad) * cos(dec * rad) * cos(lha * rad);
altitude = asin(altitude) / rad;
azimuth = sin(lha * rad) /
(cos(lha * rad) * sin(latitude * rad) -
tan(dec * rad) * cos(latitude * rad));
/* correct for quadrant */
if (lha <= 180.) {
if (azimuth > 0)
azimuth = 180. + atan(azimuth) / rad;
else
azimuth = 360. + atan(azimuth) / rad;
} else {
if (azimuth <= 0.)
azimuth = 180. + atan(azimuth) / rad;
else
azimuth = atan(azimuth) / rad;
}
printf(" %5.1f %5.1f\n", altitude, azimuth);
} else
printf("\n");
return (0);
}
double
range(val)
double val;
{
while (val < 0.0) {
val = val + 360.0;
}
if (val > 360.0) {
val = val - (to_int(val / 360.0) * 360.0);
}
return (val);
}
/*
* Helio_trans converts heliocentric ecliptical coordinates to geocentric
* right ascension and declination. It also outputs the converted value.
*/
int
helio_trans(r, b, ll, Stheta, Sr, epli, mag, sym, str)
double r, b, ll, Stheta, Sr, epli, mag;
char *sym, *str;
{
double N, D, lambda, delta, beta;
N = r * cos(b * rad) * sin((ll - Stheta) * rad);
D = r * cos(b * rad) * cos((ll - Stheta) * rad) + Sr;
lambda = atan2(N, D) / rad;
if (lambda < 0.0) {
lambda = lambda + 360.0;
}
lambda = range(lambda + Stheta);
delta = sqrt(N * N + D * D + (r * sin(b * rad)) * (r * sin(b * rad)));
beta = asin(r * sin(b * rad) / delta) / rad;
geo_trans(lambda, beta, epli, delta, mag, sym, str);
return (0);
}
/*
* Geo_trans converts geocentric ecliptical coordinates to geocentric right
* ascension and declination. It also outputs the converted value. Note
* that the output coordinates are referred to the mean equinox of date. If
* these coordinates are to be used in conjunction with a star database, use
* the epoch program to convert the planet's coordinates to the epoch for the
* star database.
*/
int
geo_trans(lambda, beta, epli, delta, mag, sym, str)
double lambda, beta, epli, delta, mag;
char *sym, *str;
{
double N, D, RA, DEC;
N = sin(lambda * rad) * cos(epli * rad)
-tan(beta * rad) * sin(epli * rad);
D = cos(lambda * rad);
RA = atan2(N, D) / rad;
DEC = asin(sin(beta * rad) * cos(epli * rad)
+ cos(beta * rad) * sin(epli * rad) * sin(lambda * rad)) / rad;
speak(RA, DEC, delta, mag, sym, str);
return (0);
}
double
htod(s)
char *s;
{
/*
* htod(str) reads floating point strings of the form {+-}hh.{mm} thus
* allowing for fractional values in base 60 (ie, degrees/minutes).
*/
double x, sign;
int full, frac;
char *t;
t = s - 1;
while (*++t) {
if (((*t < '0') || (*t > '9')) && (*t != '.') && (*t != '+') && (*t != '-'))
die("non-digit in dd.mm style numeric argument");
}
if (s == NULL)
return 0.0;
full = frac = 0;
sign = 1.0;
if (*s == '-') {
sign = -1.0;
s++;
} else if (*s == '+')
s++;
while (*s && *s != '.')
full = 10 * full + *s++ - '0';
if (*s++ == '.') {
if (*s)
frac = 10 * (*s++ - '0');
if (*s)
frac += *s++ - '0';
if (frac > 59)
frac = 59;
}
x = (double) full + ((double) frac) / 60.0;
return sign * x;
}
die(a, b)
char *a, *b;
{
fprintf(stderr, "%s: ", progname);
fprintf(stderr, a, b);
fprintf(stderr, "\n");
fflush(stderr);
exit(1);
}